#1.数据准备 

#1.1安装并加载GEOquery包（如果尚未安装）
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)
#1.2加载本地下载的系列矩阵文件
gse <- getGEO(filename = "D:/geo-database-analysis/GSE43292/GSE43292_series_matrix.txt", getGPL = FALSE)
#1.3提取表达矩阵
expr_matrix <- exprs(gse)
# 查看矩阵维度（行数=基因数，列数=样本数）
dim(expr_matrix)
# 查看前几行几列
expr_matrix[1:5, 1:3]
#1.4提取样本信息（表型数据）
pheno_data <- pData(gse)
# 查看样本信息的前几列
head(pheno_data)

#2.数据质控和检查

#2.1检查表达量分布（应该是正态或近似正态）
boxplot(expr_matrix, main = "样本表达量分布", las = 2)
#2.2检查是否有缺失值
sum(is.na(expr_matrix))
#2.3查看样本分组信息（为后续分析做准备）
colnames(pheno_data)
# 查看样本分组情况
table(pheno_data$characteristics_ch1)
#2.4对数转换（如果数据需要，芯片数据通常已经标准化）
expr_matrix_log <- log2(expr_matrix + 1)

#3.探针注释完整流程(不依赖org.Hs.eg.db包)

#3.1加载所需包
library(dplyr)
library(tidyr)
#3.2读取之前加载的GSE数据（确保确保已成功加载gse对象）
# 如果尚未加载GSE数据，请先执行：
# gse <- getGEO(filename = "你的GSE文件路径/GSE43292_series_matrix.txt", getGPL = FALSE)
# gse <- gse[[1]]
#3.3清理GPL数据，只保留有效信息
gpl_clean <- gpl %>%
  # 只保留探针ID（ID列）和基因信息列（X列）
  select(probe_id = ID, gene_info = X) %>%
  # 过滤基因信息为空或仅含空格的行
  filter(!is.na(gene_info) & gene_info != "" & trimws(gene_info) != "") %>%
  # 提取第一个参考序列ID（如NM_001004195）
  mutate(
    # 去除逗号及后面的内容，只保留第一个ID
    gene_id = sub(",.*", "", gene_info),
    # 去除可能的空格
    gene_id = trimws(gene_id)
  ) %>%
  # 过滤处理后仍为空的行
  filter(gene_id != "")
#3.4查看清洗结果
cat("有效探针-基因映射数：", nrow(gpl_clean), "\n")
head(gpl_clean[, c("probe_id", "gene_id")])
#3.5构建映射关系
probe2gene <- setNames(gpl_clean$gene_id, gpl_clean$probe_id)
#3.6筛选并转换表达矩阵
# 确保表达矩阵行名与probe2gene的命名一致（去除可能的空格）
rownames(expr_matrix) <- trimws(rownames(expr_matrix))
#3.7筛选共同探针
common_probes <- intersect(rownames(expr_matrix), names(probe2gene))
cat("表达矩阵中匹配的探针数：", length(common_probes), "\n")
#3.8生成最终表达矩阵
if (length(common_probes) > 0) {
  expr_valid <- expr_matrix[common_probes, , drop = FALSE]
  rownames(expr_valid) <- probe2gene[common_probes]
  
  # 合并重复基因ID（取平均值）
  expr_df <- as.data.frame(expr_valid)
  expr_df$gene_id <- rownames(expr_df)
  expr_final <- aggregate(. ~ gene_id, data = expr_df, FUN = mean, na.rm = TRUE)
  rownames(expr_final) <- expr_final$gene_id
  expr_final <- as.matrix(expr_final[, -1])
  
  cat("最终表达矩阵维度：", dim(expr_final), "\n")
  write.csv(expr_final, "final_matrix.csv", row.names = TRUE)
} else {
  cat("未找到匹配的探针，请检查探针ID格式是否一致！\n")
}

#4.基本的探索性分析

#4.1先从最终表达矩阵中提取样本名（定义sample_names）
sample_names <- colnames(expr_final)  # expr_final是你之前生成的表达矩阵
cat("提取的样本数量：", length(sample_names), "\n")  # 应输出64（符合实验设计）
#4.2按实验设计创建分组和患者ID
# 方案1：若样本顺序是“前32个正常组织 + 后32个斑块组织”
group <- c(rep("Normal", 32), rep("Plaque", 32))  # Normal=正常组织，Plaque=斑块组织
patient_id <- rep(paste0("Patient_", 1:32), each = 2)  # 每例患者对应2个样本
#4.3构建样本信息数据框
sample_info <- data.frame(
  Sample_Name = sample_names,    # 样本名（与表达矩阵列名完全一致）
  Patient_ID = patient_id,       # 患者ID（用于后续配对差异分析）
  Tissue_Type = factor(group),   # 核心分组（转换为因子类型，便于统计分析）
  row.names = sample_names       # 行名设为样本名，后续分析时能直接关联
)
#4.4验证样本信息是否正确
cat("\n样本信息前6行：\n")
print(head(sample_info))
cat("\n各组样本数量：\n")
print(table(sample_info$Tissue_Type))  # 应显示Normal=32，Plaque=32

#5差异分析

#5.1加载必备包
if (!require("limma")) install.packages("limma")  # 差异分析核心包
if (!require("dplyr")) install.packages("dplyr")  # 数据处理
library(limma)
library(dplyr)
#5.2构建设计矩阵（关键：纳入Patient_ID控制个体差异）
# 设计公式：表达值 ~ 0 + 组织类型 + 患者ID（消除患者间差异）
design <- model.matrix(~ 0 + Tissue_Type + Patient_ID, data = sample_info)
# 简化设计矩阵的列名（去掉前缀，更易读）
colnames(design) <- gsub("Tissue_Type", "", colnames(design))
colnames(design) <- gsub("Patient_ID", "Pat_", colnames(design))
head(design[, 1:5])  # 查看前5列（Tissue_Normal、Tissue_Plaque、前3个患者ID）
#5.3定义比较组（目标：斑块组织 vs 正常组织的差异）
contrast_matrix <- makeContrasts(
  Plaque - Normal,  # 计算“斑块组织”相对于“正常组织”的表达变化
  levels = design
)
colnames(contrast_matrix) <- "Plaque_vs_Normal"  # 给比较组命名
#5.4拟合线性模型并进行差异分析
fit <- lmFit(expr_final, design)          # 拟合线性模型
fit2 <- contrasts.fit(fit, contrast_matrix)# 应用比较组
fit2 <- eBayes(fit2)                      # 贝叶斯修正（提高统计效力）
#5.5提取差异基因结果（筛选标准：adj.P.Val < 0.05 + |logFC| > 1）
# number=Inf：提取所有基因；adjust="fdr"：用FDR校正P值
de_result <- topTable(
  fit2,
  coef = "Plaque_vs_Normal",  # 对应比较组
  number = Inf,
  adjust = "fdr",
  lfc = 1  # log2倍变化绝对值≥1（表达变化显著）
)
#5.6查看差异分析结果
cat("差异表达基因总数（adj.P.Val<0.05 & |logFC|>1）：", nrow(de_result), "\n")
if (nrow(de_result) > 0) {
  # 区分“上调基因”（斑块中表达更高）和“下调基因”（斑块中表达更低）
  up_genes <- de_result %>% filter(logFC > 0) %>% rownames()
  down_genes <- de_result %>% filter(logFC < 0) %>% rownames()
  cat("上调基因数：", length(up_genes), " 下调基因数：", length(down_genes), "\n")
  # 显示前10个差异基因
  print(head(de_result, 10))
} else {
  cat("未筛选到显著差异基因（可放宽标准，如adj.P.Val<0.1或|logFC|>0.5）\n")
  # 若无显著基因，放宽标准重新提取
  de_result <- topTable(fit2, coef = "Plaque_vs_Normal", number = Inf, adjust = "fdr", lfc = 0.5)
}

#5.7保存差异分析结果（方便后续查看和论文使用）
write.csv(de_result, "差异表达基因结果.csv", row.names = TRUE, fileEncoding = "UTF-8")

#6.可视化分析

#6.1火山图
# 6.1.1基于差异分析结果创建volcano_data（包含基因ID）
volcano_data <- as.data.frame(de_result)  # 将差异结果转换为数据框
volcano_data$Gene_ID <- rownames(de_result)  # 手动添加基因ID列（原行名）
# 6.1.2确认volcano_data已创建
cat("volcano_data是否存在：", exists("volcano_data"), "\n")  # 应显示TRUE
head(volcano_data[, c("Gene_ID", "logFC", "adj.P.Val")])  # 查看关键列
# 6.1.3添加分组标签（完全不依赖其他包）
volcano_data$Group <- ifelse(
  volcano_data$adj.P.Val < 0.05 & volcano_data$logFC > 1,
  "显著上调",
  ifelse(
    volcano_data$adj.P.Val < 0.05 & volcano_data$logFC < -1,
    "显著下调",
    "非差异"
  )
)

#6.1.4查看结果
head(volcano_data[, c("Gene_ID", "logFC", "adj.P.Val", "Group")])
#6.1.5绘制火山图
ggplot(volcano_data, aes(x = logFC, y = -log10(adj.P.Val), color = Group)) +
  geom_point(size = 2, alpha = 0.7) +  # 散点图（alpha=透明度）
  #设置颜色（上调红、下调蓝、非差异灰）
  scale_color_manual(values = c("显著上调" = "#E74C3C", "显著下调" = "#3498DB", "非差异" = "#95A5A6")) +
  # 添加阈值线（虚线：P=0.05和logFC=±1）
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black", size = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black", size = 0.5) +
  # 标题和标签（适配论文格式）
  labs(
    x = "log2(表达倍数变化) [logFC]",
    y = "-log10(校正后P值) [-log10(adj.P.Val)]",
    title = "颈动脉斑块组织 vs 正常组织 差异表达基因火山图",
    color = "基因分类"
  ) +
 
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.position = "bottom"
  )
#6.1.6保存火山图
ggsave("火山图.png", width = 10, height = 8, dpi = 300)

#6.2热图

#6.2.1下载pheatmap包
if (!require("pheatmap")) install.packages("pheatmap")
library(pheatmap)
#6.2.2准备数据：选择前20个差异最显著的基因（不足则取全部）
n <- min(20, nrow(de_result))  # 最多显示20个基因
if (n == 0) stop("无差异基因可用于绘制热图")
#从表达矩阵中提取这些基因的表达值
heatmap_genes <- rownames(de_result)[1:n]
heatmap_data <- expr_final[heatmap_genes, ]
# 数据标准化（按基因标准化，使不同基因的表达趋势可比）
heatmap_data_scaled <- t(scale(t(heatmap_data)))
# 准备样本分组注释（用于热图上方标记）
sample_annot <- data.frame(
  Tissue_Type = ifelse(
    sample_info$Tissue_Type == "Normal", 
    "正常组织", 
    "斑块组织"
  ),
  row.names = rownames(sample_info)
)

#6.2.3绘制热图
pheatmap(
  mat = heatmap_data_scaled,
  annotation_col = sample_annot,  # 样本分组注释
  annotation_colors = list(
    Tissue_Type = c("正常组织" = "#4285F4", "斑块组织" = "#EA4335")
  ),
  show_rownames = TRUE,  # 显示基因ID
  show_colnames = FALSE, # 隐藏样本名（避免拥挤）
  treeheight_row = 15,   # 基因聚类树高度
  treeheight_col = 15,   # 样本聚类树高度
  color = colorRampPalette(c("#3366CC", "white", "#DC3912"))(100),  # 蓝-白-红配色
  main = paste("Top", n, "差异基因表达热图")  # 标题
)

#6.2.4保存热图
dev.copy(png, "差异基因热图.png", width = 1000, height = 800, res = 300)
dev.off()

#6.3差异基因数量柱状图

#6.3.1检查差异分析结果是否存在
  if (!exists("de_result") || nrow(de_result) == 0) {
    stop("❌ 差异分析结果无效，无法绘制柱状图")
  }
  
  # 重新统计差异基因（使用基础R函数，不依赖dplyr）
  # 显著上调：adj.P.Val < 0.05 且 logFC > 1
  up_count <- sum(de_result$adj.P.Val < 0.05 & de_result$logFC > 1, na.rm = TRUE)
  # 显著下调：adj.P.Val < 0.05 且 logFC < -1
  down_count <- sum(de_result$adj.P.Val < 0.05 & de_result$logFC < -1, na.rm = TRUE)
  
  # 处理无显著差异基因的情况
  if (up_count == 0 & down_count == 0) {
    cat("⚠️ 未检测到显著差异基因，放宽筛选标准...\n")
    up_count <- sum(de_result$adj.P.Val < 0.1 & de_result$logFC > 0.5, na.rm = TRUE)
    down_count <- sum(de_result$adj.P.Val < 0.1 & de_result$logFC < -0.5, na.rm = TRUE)
  }
  
  cat("差异基因统计：上调", up_count, "个，下调", down_count, "个\n")
#6.3.2柱状图数据
  bar_data <- data.frame(
    类型 = c("显著上调", "显著下调"),
    数量 = c(up_count, down_count)
  )
  
  # 设置保存路径
  bar_path <- "D:/geo-database-analysis/GSE43292/image/差异基因数量柱状图.png"
  
  # 显式创建图片设备
  png(
    filename = bar_path,
    width = 1000,
    height = 800,
    res = 300,
    family = "SimHei"  # 支持中文
  )
  
#6.3.3绘制柱状图
  bar_heights <- bar_data$数量
  bar_names <- bar_data$类型
  
  barplot(
    height = bar_heights,
    names.arg = bar_names,
    col = c("#EA4335", "#4285F4"),
    main = "差异表达基因数量统计",
    xlab = "基因表达变化类型",
    ylab = "基因数量",
    ylim = c(0, max(bar_heights) * 1.2),  # 预留顶部空间
    cex.main = 1.2,
    cex.lab = 1.1,
    cex.names = 1.1
  )
  
  #在柱子上添加数字标签
  text(
    x = 1:length(bar_heights),
    y = bar_heights + max(bar_heights) * 0.05,
    labels = as.character(bar_heights),
    cex = 1.2,
    font = 2
  )
  #关闭图片设备（关键步骤）
  dev.off()

#6.3.4验证是否生成文件
  if (file.exists(bar_path)) {
    cat("✅ 差异基因数量柱状图已保存至：", bar_path, "\n")
  } else {
    warning("⚠️ 差异基因数量柱状图保存失败，请检查路径是否存在")
  }
  
  
  
  


  
 